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The dynamic stability of the Boolean networks representing a model for the gene transcriptional regulation 
(Kauffman model) is studied by calculating analytically and numerically the Hamming distance between two 
evolving configurations. This turns out to behave in a universal way close to the phase boundary only for 
in-degree distributions with a finite second moment. In-degree distributions of the form Pd(k) ~ k~"* with 
2 < y < 3, thus having a diverging second moment, lead to a slower increase of the Hamming distance when 
moving towards the unstable phase and to a broadening of the phase boundary for finite N with decreasing y. We 
conclude that the heterogeneous regulatory network connectivity facilitates the balancing between robustness 
and evolvability in living organisms. 

PACS numbers: 89.75.Hc,64.60.Cn,05.65.+b,02.50.-r 
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I. INTRODUCTION 

Complete genome sequencing and the analysis of the bind- 
ing of transcriptional regulators to specific promoter se- 
quences have uncovered the global organization of the gene 
transcriptional regulatory network in well-studied organisms 
such like Escherichia coli 01 and yeast Saccharomyces cere- 
visiae [2]. The gene network describes a directed relationship 
- regulation - between different genes, and its architecture is 
characterized by broad connectivity distributions U |H S 0], 
over-representation of selected motifs JH, and so on. These 
features are rarely found in random networks, probably being 
the consequence of evolutionary selection. Therefore illumi- 
nating the functional characteristics associated with those dis- 
covered structural features can help trace back their origins. 
In this work, we show heterogeneous connectivity can facil- 
itate the balancing between dynamical stability and instabil- 
ity. Both robustness and evolvability are essential for living 
organisms, which achieve their specific phenotype by then- 
gene expression program [6]. Thus the transcriptional reg- 
ulatory network should be organized in a way that supports 
the coexistence of these apparently contradictory properties 
and from this perspective, it has been proposed that the gene 
network should be at the boundary between stable and unsta- 
ble phases, called the edge of chaos 01. The question then 
arises: What are the characteristics of the network architec- 
ture that can support the requirement to be located at the edge 
of chaos? A simple model incorporating recently available 

information turns out to be useful to answer this question. 

r— i 

The Kauffman model QVTJ was used in the past to study 
the gene network dynamics which is far from completely 
known because of its complexity. In this model, each node 
has a Boolean variable, 1 or 0, the discretized expression 
level, evolving regulated by other K nodes according to the 
quenched rules that are randomly distributed with a param- 
eter p. In spite of these simplifications involved, the model 
revealed detailed relations between the dynamical stability 
against perturbations and the network architecture Si- 
Moreover, distinct attractors in the configuration space are 



considered as corresponding to different cell types in a given 
organism and thus its scaling with the number of genes 
(nodes) across different organisms has been of great inter- 
est 0, § E E E E Q. Empirically, the number of 
cell types scales as the square-root of the number of genes, 
and the same scaling relation was believed to hold between 
the number of attractors and the number of nodes in the 
Kauffman model at the critical point p c (K), supporting the 
hypothesis that living organisms should be between order 
and chaos [7]. Recently, however, it was found that under- 
sampling effects may hamper numerical enumeration of dis- 
tinct attractors ifTol [Till and further investigations demon- 
strated that the total number of attractors grows faster than any 
power law with the system size Il2l 1 1 3fl . On the other hand, it 
was also reported that attractors stable against deviation from 
synchronous update show a sub-linear scaling behavior fl4ll . 

Recent investiations of real gene networks suggest general- 
ization of the original Kauffman model. First, the distribution 
of the regulating rules is structured showing a bias towards the 
canalyzing functions |15Ulql . Second, the number of links or 
degree is not constant but different from node to node, result- 
ing in broad degree distributions lfl7ll . In the gene regulatory 
networks of E. coli HUE and yeast USEE, El, the 
distributions of out-degree (number of target genes for each 
regulator) and in-degree (number of regulators for each tar- 
get gene) were not delta-functions but shown to take power- 
law or exponential-decaying form, respectively, although true 
asymptotic behaviors were hard to discern due to finite size 
effects. While the effects of the structured distribution of reg- 
ulating rules have been intensively studied lfl6[ E. El, it 
remains to show how the heterogeneous connectivity affects 
the dynamical stability |EE- 

We consider the Kauffman model on directed networks 
with general in- and out-degree distributions and compare two 
evolving dynamical configurations by computing their Ham- 
ming distance, to determine whether a given network is dy- 
namically stable (zero distance) or unstable (non-zero dis- 
tance) against perturbations. The critical point of the Boolean 
networks with power-law degree distributions was studied re- 
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cently ll24ll . In the present work, we show quantitatively how 
the Hamming distance behaves near the critical point, which 
will provide a deeper understanding of the critical phenom- 
ena of Boolean networks with heterogeneous connectivity pat- 
terns and insights into the interplay of structure and dynam- 
ics in living organisms. The Hamming distance for infinite 
system size (thermodynamic limit) can be computed by the 
method presented in Ref. [18] and we here present a detailed 
description of the method along with a discussion on the ef- 
fects of correlation between in- and out-degree of the same 
node. Then, more importantly, we extend the method to derive 
the Hamming distance for finite system size, which enables us 
to check the analytic predictions with numerical simulation re- 
sults. Our main result is that for in-degree distributions with 
a diverging second moment the Hamming distance increases 
very slowly when moving from the phase boundary towards 
the unstable phase and the width of the boundary in finite-size 
systems is very broad. This indicates that strongly heteroge- 
neous genetic networks have a large capacity to stay at the 
edge of chaos when their structural and functional organiza- 
tion is subject to variation. 

The paper is organized as follows. We introduce the Kauff- 
man model for Boolean networks in Sec. [II] In Sec. [TTTJ the 
annealed approximation is described and used to compute 
the Hamming distance, which reveals different phases of the 
Boolean networks. The finite-size effects on the critical phe- 
nomena of the Boolean networks are derived using the an- 
nealed approximation in Sec. [IV] Finally, the results are sum- 
marized and discussed in SecfVl 



II. MODEL AND HAMMING DISTANCE 

In the Kauffman model, the dynamical configuration of N 
Boolean variables at time f, E(f) = {o,(f)|/ = 1,2,... ,N}, is 
updated in parallel as 

a,(r + l) =/,•(£,•(')), (1) 

where a, for each i takes 1 or and E,-(?) = 
{oi l (t),Oi 2 (t),...,Oi,(t)} denotes the configuration at 
time t of the ki regulators %s = {h,i2, ■ ■ ■ ,hj}, of the 
node i. The functional dependencies between nodes via 
{fi(T<i) | i = 1,2,..., N} constitute a directed network in which 
two nodes i and j are connected with a directed edge if 
j 6 Hi, where is an outgoing edge of node j and an in- 
coming edge of node i. The quenched, i.e., time-independent, 
regulating rules are random Boolean functions, i.e., they are 
chosen randomly such that /;(£;) for a given E, is 1 with 
probability < p < 1 and with probability 1 — p. The 
parameter p deviating from 1/2 indicates an asymmetry 
between expressed (1) and non-expressed (0) state of a gene. 

We focus on the following question: If one starts at time 
t = with two randomly chosen configurations, E and E with 
<Sj 7^ Oj for all j, that is, all node states perturbed (altered), 
how many nodes remain perturbed at time t > 0? The fraction 
of these perturbed nodes or the Hamming distance between E 



and E at time t is defined as 

iv i=l 

with 8 a j, being 1 for a = b and otherwise. The Hamming 
distance may vary between and 1 depending on the dynamic 
asymmetry parameter p. We will see in the next section that 
the value of the Hamming distance in the stationary state may 
display a transition from zero to a non-zero value as the net- 
work architecture and the parameter p are varied. 



III. ANNEALED APPROXIMATION AND PHASE 
TRANSITION OF BOOLEAN NETWORKS 

In this section, we investigate the phase diagram of the 
Kauffman Boolean network defined in the previous section 
by computing analytically and numerically the Hamming dis- 
tance for infinite system size. This allows us to understand dif- 
ferent phases of the Boolean networks determined by network 
structure and the parameter of dynamic asymmetry. Some of 
the results presented in this section are also found in Ref. Hill . 



A. Annealed approximation 

A recursion relation for the Hamming distance between 
consecutive time steps is obtained by the "annealed" approxi- 
mation [8]. While the regulation rule /,■ and the regulators % 
are fixed for each node i in the original model, one assigns 
them randomly to every node at every time step, keeping the 
in-degrees and the out-degrees, in the annealed approxima- 
tion. Then the evolution of the Hamming distance H^g^t) 
for the nodes with in-degree k and out-degree q is given 
by H Kq {t + 1) = X[l - (1 - L k ', q >q'P<i(k',q')H k , tg ,(t)/{q)) k }, 
where X = 2p(l — p), Pd(k',q') is the joint distribution of 
k' and q' , and (q) = Y.k,q c lPd{k 1 q)- The correlation between 
the degrees of neighboring nodes, {k,q} and {k',q'} is ig- 
nored in this formalism but will be discussed in Sec. [V] 
The parameter X ranging from to 1/2 is the probabil- 
ity that yields different outputs for E; and E, different 
and the term within the brackets represents the probabil- 
ity of the latter. Note that the degree distribution for the 
regulators is weighted by their out-degrees. If we intro- 
duce H(t) = Y l k,tj[qPd(k,q) I (q)]Hk,q(t), it is obtained self- 
consistently and in turn H(t) — Y l k,qPd(k,q)Hi ( ,q(t) is com- 
puted as follows [18]: 

H(t + l) = 

H{t + \) = X^Pdim- (1-^(0)*]. ( 3 ) 

k 

where P c i(k) — Y.qPd(k,q) is the in-degree distribution. In the 
original Kauffman model where the in-degree is fixed to k = 
K, Eq. ©reduces to H(t+1) =X[l -(1 -H(t)) K ] d. 
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FIG. 1 : (Color online) Hamming distance H for the Kauffman model 
as a function of X. The points are the simulation results for model net- 
works H3] withW= 10 4 and (k) = 4 that have P d (k,q) = gy{k)g n (q), 
where g fl (x) ~ x~ a for a finite and goo(x) = (k) x e~W /x\. The Ham- 
ming distance H was first averaged over time between t = 100 and 
/ = 200, and then averaged over 1000 different realizations of net- 
works and regulating rules for which H ^ 0. The lines are the nu- 
merical solutions to Eq. ([3} with H(t) = H(t) = H for all t and the 
same degree distributions as in the simulation used. (Inset) Plots of 
InH versus ln(X — X c ) with X c = 0.25. The slopes agree with Eq. l[8}. 



B. Ordered and chaotic phases 

The limiting value H = lim r ^oo H(t) can characterize the 
system's response to dynamical perturbations. Replacing H{t) 
and H(t) with H and H, respectively, and expanding the first 
line in Eq. (0 for small H, one finds that H = and H = 
for X < X c and H > and H > for X > X c ■, where the critical 
point X c depends on the network topology via 



Xr = K 



-1 



with K : 



k.q 



kqP d (k,q) 
(q) ^ 



(4) 



That is, the system is in the ordered phase, any perturbation 
making no effect on the system eventually, when X <X C . On 
the other hand, the system does not remain in the same sta- 
tionary state but shifts to another stationary state triggered by 
a perturbation when X > X c . 

The quantity K may be considered as the average in-degree 
of regulators weighted by their out-degrees. Introducing the 
conditional average of out-degree, = Y,qlPd(k,q)/Pd{k), 
we can rewrite the quantity K as K — £ t k(qk / ' (q))P d (k) and 
also the first relation in Eq. (0 as 



H(t + l)=X 1 £^-P d (k)[l-(l-H(t)) k ]. 



(5) 



If qt is independent of k or (more strongly) the in-degree and 
the out-degree of a node is not correlated statistically, it fol- 
lows that qk = (q) for all k, K reduces to the conventional av- 
erage in-degree (k) = Y<kkPd(k), and H(t) — H(t). The anal- 
yses of the transcriptional regulatory networks of E. coli and 
yeast show no significant variation of q^ with k [25] and so 
we will assume in the following that q^ = (q) for all k. In 



Sec. lIVCl we will discuss how our results for the critical phe- 
nomena would be changed by the ^-dependence of q^. Under 
this assumption (q^ = (q)), the Hamming distance H{t) de- 
pends only on the in-degree distribution P,i (k) and the dynam- 
ics parameter A,. 



It has been shown that the scaling behavior of the average 
number of attractors with the system size remains to be the 
same for different out-degree distributions such as uniform, 
exponential, and power-law one 12311 . Our analysis based on 
the annealed approximation suggests further the irrelevance of 
the out-degree distribution to the Hamming distance. To con- 
firm this as well as check the validity of the annealed approxi- 
mation or Eq. (0, we performed simulations of the Kauffman 
model defined in Season an ensemble of model networks 
constructed as follows 112611 : i) Each of nodes has two in- 
dices ;'i n and / out , which run from 1 to jV respectively. The 
two indices are given independently to each node, ii) Choose 
a node A with index / ou t with probability j^ u 't 0Ut /E/ J™ *- iii) 
Choose a node B indexed i m with probability / in am / X; 7 ain • 
iv) Assign a link from the node A to B unless they are con- 
nected, v) Repeat ii) and iii) until the total number of links is 
L. The generated networks have nodes, L links, and degree 
distribution given by P c j(k,q) — P c j(k)P c i(q), where in-degree 
distribution takes the form P c i {k) ~ k~1 and the out-degree dis- 
tribution takes the form P c i{q) ~ q~^ with y = 1 + 1 /a m and 
T| = 1 + 1/ a out . It is then obvious that q^ = (q) for all k. When 
0Cj n = and thus all the nodes can have an incoming link with 
equal probability, the in-degree distribution becomes a Pois- 
son one, P d (k) = (k) k e-^/kl. Thus the de gree distribution 
may take power-law or Poissonian form depending on the val- 
ues of 0Cj n and a ou t corresponding to scale-free (SF) networks 
or completely random networks, respectively. The simulation 
results (data points) shown in Fig. Q] are compared with the 
numerical solutions (lines) to Eq. (0, the annealed approxi- 
mation, which show a good agreement and support the valid- 
ity of the annealed approximation. Also it is shown that the 
Hamming distance is the same for different out-degree distri- 
butions. 



The implication of Eq. (0 for SF networks has been dis- 
cussed in Ref. 3, where P d (k) = Ar7t;(Y) for k = 1,2,... 
with f^(x) the Rieman-zeta function. Based on the result 
that (k) = C(Y- 1)/C(Y) < 2 for y > 2.47875..., it was 
claimed 12411 that the abundance of SF networks with 2 < y < 
2.5 in nature and society can be attributed to the presence of 
both phases, stable and unstable, only in such networks. How- 
ever, the values of (k) and y do not show such strong correla- 
tion in real networks. For instance, the average degree (k) 
ranges from 2.57 (Internet router network) to 28.78 (movie 
actor network) although the degree exponent y lies between 2 
and 3 II 1711 , which is possible due to the power-law behavior 
observed only asymptotically. We will show in the next sec- 
tion that root for the dynamical advantage of SF network lies 
elsewhere. 
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C. Critical exponents 

Here we address the behavior of the Hamming distance 
around X c for infinite system size. In that regime of X, the 
Hamming distance is very small and its increase with X — X c: 
can be characterized by a scaling exponent. This critical be- 
havior of the Hamming distance is of interest to us because 
it shows how the network topology is related to the system's 
dynamic response. 

When the moments {k") = Y,k p d(k)k n for all n > are fi- 
nite, Eq. (0 can be written as 



H 



'.X^^—i^H". 



(6) 



Keeping the leading terms, we find that H ~ (X/X C )H — 
X(k 2 )H 2 /2, which gives 

H ~A 

with A = X/X c — 1 for < A <C 1. This result can be repre- 
sented as H ~ A^ with (3 = 1, where we introduced the critical 
exponent p. 

On the other hand, if Pd(k) ~ ck~^ for k ^ 1 with c a con- 
stant, {k n } diverges as ck^l ^ / (n + 1 — y) for n > fy] — 1 with 
^max the largest in-degree and \x] denoting the smallest inte- 
ger not smaller than x. Applying the relation Y,k>k max Pd(k) ~ 
\/N from the extreme value statistics 02711 . one can see that 
£ max scales as ATVu-l), The diverging terms in Eq. (|6]l 
have alternative signs and lead to non-analytic terms in H 
as described below [28]. For small H, Eq. (f3j) reads as 
H ~ XY,kPd(k)[l — e~ kH ] and recalling the power-law form 
of Pd(k), Pd(k) ~ ck^, we can utilize the fact that the Mellin 
transform of a function F(y,H) = Y*k=i k~^e~ kH is given by 
^(Y,*) = So F {l,l£)H s - x dH = T{s%{s+i) with T{x) the 
Gamma function Il28ll . The inverse transform of f (y,s) 
then is represented in terms of the poles of the Rieman-zeta 
function and the Gamma function, which gives F(y,H) = 

Th,s)H-*ds = T{l-i)W- l + ^: =0 {-lYln\Uy- 
n)H" [28]. Therefore we find that, for P d {k) ~ ckr^ |29||, 



rvi-2 



(-1) 



B+l 



■ (k n )H" -XcT(l- y)H^ 1 + o (// M - 1 ) . 

(7) 

If y > 3, the H 2 term is the next leading term in the right-hand- 
side of Eq. (O and then the critical behavior is given by H ~ A 
as in the case of all (k") finite. On the other hand, if 2 < y < 3, 
the W~ 1 term becomes the next leading term and we find that 
H ~ (X/X C )H - XcF(l - f)W-\ Therefore for X > X c , 

H~A l /(y- 2 \ 

In summary, we can list the values of the critical exponent P 
varying with the in-degree exponent y as lfl8ll 

r _ / 1 (Y>3), 

p -\ l/(y-2) (2<y<3), W 

which is confirmed numerically [See the inset of Fig. [Q- 




FIG. 2: (Color online) Phase diagram of the Kauffman model for (a) 
a Poisson in-degree distribution and (b) a power-law one with the 
exponent y = 2.5, both in case of uncorrelated in- and out-degree. 
The color encodes the Hamming distance H obtained by numerically 
solving Eq. ifj} under setting H(t) = H(t) = H. 
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FIG. 3: Finite-size scaling behavior in the critical regime, (a) Data 
collapse for different system sizes with y — » °°. The inset shows that 
Xc ~ 0.2505(5). (b) Data collapse with y = 2.5 and X c ~ 0.287(1). 



The critical exponent p varying with the in-degree distri- 
bution as in Eq. © is illuminating how the network topology 
affects the system's response to perturbation. Figure |2] shows 
phase diagrams of the Kauffman model on model networks 
with a Poisson in-degree distribution and with a power-law in- 
degree distribution with the exponent y = 2.5, in which color 
represents the value of the Hamming distance. As shown in 
the figure, the SF networks with 2 < y < 3 and thus larger 
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values of (3 keep the Hamming distance non-zero but small 
in a much larger region in the (A, (k)) plane than those with 
y > 3. Structural and functional organization of cellular net- 
works, parameterized here by (k) and A (/?), respectively, may 
be subject to unexpected changes. Our finding suggests that 
the systems with strong heterogeneous connectivity patterns 
can maintain their dynamic criticality robustly, and further, 
votes for the hypothesis that living organism's machinery lies 
at the edge of chaos. 



IV. BOOLEAN NETWORKS OF FINITE SIZE 

In the thermodynamic limit N — > °°, the critical point does 
not exhibit any dependence on the network topology. How- 
ever, for finite N, the critical point itself develops its depen- 
dence on the network topology: It is no more a point but has 
a non-zero width depending on N. Adopting the finite-size 
scaling ansatz j30ll 



(9) 



with the scaling function *¥(x) — * const, for x -C 1 and 
^(x) — > x$ for x> 1, one can see that H ~ in the crit- 

ical regime \AN l ^\ < l, In the A, axis, this critical regime 
is N~ 1 ^ wide for N finite and shrinks to zero in the thermo- 
dynamic limit. Therefore the scaling exponent /j describes 
the width of the critical regime for finite-size systems. In the 
critical regime, the cluster of perturbed nodes, explained be- 
low, exhibits scale invariance characterized by a power-law 
distribution of its size, which is connected to the behavior 
H ~ N~P/ M . We show in the next that the asymptotic behavior 
of the cluster size distribution can be derived using Eqs. © 
and @, which allows us to obtain the scaling exponent /j and 
to check Eq. (|9j. 



A. Evolution of perturbed-node clusters 

The parameter A, denotes the probability that a node be- 
comes perturbed (o,(f + 1) ^ d;) once the configuration of 
its neighbors are perturbed (£,■ ^ £,-). When A is zero, the 
Hamming distance becomes zero immediately even though all 
nodes were perturbed initially. As X increases from 0, clusters 
appear, consisting of perturbed nodes that are connected by 
active edges. An edge from node A to B is inactive if the per- 
turbation at node A cannot bring any difference to the dynami- 
cal state of node B, and active otherwise. While a node's state 
can be totally irrelevant to its neighbor connected by an inac- 
tive edge, perturbation at one node can propagate to its neigh- 
bors through active edges. The perturbed-node clusters evolve 
with time, decaying or growing. A perturbed node may be- 
come normal (non-perturbed) due to its regulators becoming 
normal at a certain time step, which leads to the decay of the 
cluster it belonged to. On the other hand, normal nodes may 
become perturbed due to its regulators becoming perturbed, 
which leads to the growth of a cluster. 



When the parameter X reaches the critical regime or goes 
beyond it (X > X c ), a giant cluster of perturbed nodes ap- 
pears and contributes to the Hamming distance in the station- 
ary state. There may exist many smaller clusters, but they are 
hard to survive eventually. A small cluster has a relatively 
small number of perturbed nodes, and they are surrounded by 
many normal nodes. Therefore the perturbed nodes in smaller 
clusters have higher chance of becoming normal than those in 
larger clusters, which leads to the higher chance of shrinking 
and decaying of smaller clusters. The perturbed nodes in the 
giant cluster, on the other hand, have more perturbed nodes as 
regulators and thus the giant cluster has a higher chance to sur- 
vive; the probability becomes nonzero when X > X c . There- 
fore the Hamming distance H can be approximated using the 
size of the largest cluster S via H ~ S/N 113 111 . We will use 
this relation in the below to derive the asymptotic behavior 
of the size distribution of the clusters from the self-consistent 
equations for H, Eqs. (|6]l and (0. 



B. Cluster-size distribution in the annealed approximation 

Let us denote the probability that a node belongs to a sizes 
cluster (of perturbed nodes) by P(s) and consider its generat- 
ing function CO = fP(z) = Y. s P{ s )z s - The Hamming distance is 
related to the generating function as 



H 



S 
N 



s<S 



: Um[l-*(4) 



(10) 



where z* N = e~ l l s and S satisfies S2 <C S <C S with 52 the sec- 
ond largest cluster size. This relationship, combined with 
Eqs. © and (0, gives us the functional form of the in- 
verse function z = fP (co). Let's expand the inverse function 
around co=lasz=l — Y.n>\ ^«(1 — <*>)"■ Then one can see 
that this equation should reduce to Eqs. © or (0, depend- 
ing on the in-degree distribution, with z and 1 — CO replaced 
by Zx and H, respectively, in the thermodynamic limit. Note 
that z% — > 1 in the thermodynamic limit. Therefore the inverse 
function z = fP ~ (co) is expanded around CO = 1 as follows: 



1 



( 1 - A, Ac) ( 1 - co) + X(k 2 ) ( 1 - co) 2 /2 + ■ 



(11) 



for y > 3 and 

l-z-(l-A/A,)(l-co) + Acr(l-Y)(l-co) Y - 1 + --- (12) 

for 2 < y < 3, where c is the coefficient appearing in the 
asymptotic behavior of the in-degree distribution Pd(k) ~ 

ck-y. 

The functional behavior of ?P (z) around z = 1 is then de- 
rived from Eqs. ( fTTT i and ( TT2l . At the critical point (A = A c ), it 
exhibits a singularity depending on the in-degree exponent: 

1 Wz)~I (1 " Z)1/2 (y>3) (13) 
1 * {Z) \ (l-z)V(Y-i) (2<Y <3) (13) 

The asymptotic behavior of the cluster-size distribution P(s) 
can be obtained from the functional form of T (z) since P(s) — 
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(l/s\)(d s /df)¥ (z)\ z =o. In particular, when 1 -${z) ~ (1 - 
z) x_1 with x a non-integer, the cluster-size distribution takes a 
power-law form P(s) ~ s~ z . The asymptotic behavior of P(s) 
is thus distinguished between y > 3 and 2 < y < 3 as follows: 



,- 3 /2 (y >3), 

(2<y<3) 



(14) 



Scale invariance is typical of the system at the critical point 
and the power-law form of the cluster-size distribution is an 
example I3TI1 . Our finding shows that the power-law exponent 
may vary with the degree exponent. In the subcritical {X < 
X c ) or supercritical (X > X c ) regime, the linear term in 1 — 
CO in Eqs. ( fTTb and ( fT2l is dominant around co — 1 and the 
cluster-size distribution takes an exponential-decaying form 
like P(s) ~ e~ s ^ 



C. Scaling exponents 

Once the asymptotic behavior of P(s) is known, the scal- 
ing property of the largest cluster size may be derived. Using 
Eq. (O in the relation I S>S P(» ~ S/N, a part of Eq. JlOb . 
one can see that the size of the largest cluster S scales with 
the system size N as S ~ N 2 / 3 for y > 3 and S ~ M? -1 ^ 
for 2 < y < 3. It is known that the number of nonfrozen 
nodes in critical Boolean networks with any fixed number 
of inputs 11321 |33[] or fast-decaying in-degree distribution also 
scales as N 2 / 3 113411 . A node i is called nonfrozen if its state 
(a,) is not fixed at either or 1 in the stationary state and thus 
may be perturbed (but not necessarily). The set of perturbed 
nodes is thus a subset of the set of nonfrozen nodes and it is 
noteworthy that their sizes display the same scaling behavior. 

The Hamming distance, H ~ S/N, in the critical regime is 
then given by 



H 



A^ 1 / 3 (Y>3), 
N- 1 ^ (2<y<3). 



(15) 



Since H ~ N~^^ in the critical regime from the finite-size 
scaling ansatz in Eq. (O, we find that the scaling exponent /j 
is given by 



^-{y/( Y - 



(Y>3), 
(2<y<3) 



(16) 



To check numerically the derived finite-size scaling be- 
havior of the Hamming distance, we performed simulations 
of the Kauffman model on the model networks described in 
Sec. IIIIBl varying the system size. First, we plot the data of 
as a function of X for the same average degree ( (k) = 4) 
and different system sizes (N = 4000,8000,16000,32000). 
They cross at one point, which determines the critical point 
X c according to Eq. (O (See the insets of Fig . [3j . Using these 
values of X c , we plotted the same data of versus AN 

in Fig. [3] The collapse of the data from different system sizes, 
while a slight deviation is seen in case of SF networks, pre- 
sumably due to strong finite-size effects, supports the scaling 



behavior of the Hamming distance in Eq. (O with Eqs. (0 and 
( [T6| > used. 

The finite-size scaling behavior in Eq. (O has been identi- 
fied also in a wide range of dynamical systems on complex 
networks. In particular, the formation of a giant cluster in SF 
networks evolving by adding links has been analyzed through 
its exact mapping to the q = 1 Potts model and found to be 
characterized by (p = 1 , fj = 3) for the degree exponent y > 4 
and (P = l/(y-3), fi = (y- l)/(y-3)) for 3 < y < 4 fU. 
These exponents are very similar to Eqs. (O and ( fTBT ). Also 
in the Ising model 11351 l36l 1371 13811 a nd the Kuramoto model 
for synchronization phenomena 11391 14011 . the exponents are 
given by(p = l/2, /j = 2) for the degree exponent y > 5 and 
(P=l/(y-3),Ai=(y-l)/(Y-3))for3<y<5. Such sim- 
ilar dependence of the scaling exponents on the degree expo- 
nent suggests a common framework to understand the critical 
phenomena on complex networks I37ll4lll42ll . 

We note that while the scaling plot gives X c ~ 0.251 for the 
completely random networks (y — > 00) as predicted by Eq. (|4j, 
the SF networks with y = 2.5 have X c ~ 0.289, deviating from 
the predicted value 0.25. This deviation seems to be rooted 
in the use of the annealed approximation for the Hamming 
distance. It has been reported that the analytic prediction 
of the critical point in the framework of the mean-field the- 
ory deviates slightly from the numerical analysis in the Ising 
model [38] and in the Kuramoto model for synchronization 
phenomena |j40] . An improvement can be made by consid- 
ering the Cayley tree with a give n degree distribution as the 
underlying network topology I14lll . 

Our results show that the width of the critical regime W ~ 
N^if in the A, axis increases as the in-degree exponent y de- 
creases below 3 while its scaling behavior remains the same 
for all y > 3. Since the number of genes in most organisms is 
not infinite but of order 10 5 at most, the width of the critical 
regime may be between o(10~ 2 ) and 0(1), depending on the 
in-degree exponent. Such a broad critical regime for small val- 
ues of y should help living organisms to remain in the critical 
regime and in turn, to balance between robustness and evolv- 
ability. Individual dynamical responses depend on the prop- 
erties of the perturbed elements, i.e,, on their connectivities 
and regulating rules, which leads to perturbation propagation 
on various scales in heterogeneous networks 82411 . Here we 
have analyzed the whole ensemble of such differentiated dy- 
namical responses in heterogeneous networks and found that 
it can remain critical more easily with the help of extremely 
heterogeneous connectivity patterns. 

The effects of correlation between in- and out-degree of the 
same node, which we have ignored so far, can be addressed 
in the results we obtained. In presence of the in- and out- 
degree correlation, (q^/ (q))P l j(k) should be considered in- 
stead of P c i(k) to compute the Hamming distance, as seen in 
Eq. (0. If it holds that qt ~ AT 9 for large k, we have to con- 
sider the effective in-degree exponent, y e ff = y+0, in place 
of y for power-law in-degree distributions in all the results we 
obtained, including the scaling exponents in Eqs. (JHJ and ( TT~6b . 
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V. SUMMARY AND DISCUSSION 

In summary, we investigated the phase transition between 
the stable (ordered) and unstable (chaotic) phase in the 
Boolean dynamical network. Heterogeneous connectivities 
are found to broaden substantially the small Hamming dis- 
tance region close to the phase boundary by suppressing the 
perturbation propagation in the unstable phase. Furthermore 
the transition region for finite system sizes turns out to be 
much wider than in homogeneous networks. Such a robust 
pseudo-criticality is expected to be also present in transcrip- 
tional regulatory networks and also in other biological net- 
works such as neural networks [43], which can be a source for 
stability and evolvability coexisting in living organisms. 

Our results suggest that the heterogeneous connectivity pat- 
terns of many biological networks have been selected in the 
course of evolution in part to serve for achieving stability 
and evolvability simultaneously. Therefore it would be de- 
sirable to propose a model in which heterogeneous connec- 
tivity patterns emerge driven by the evolutionary pressure to- 
wards a broad edge of chaos. There are indeed models for 
co-evolution of network structure and dynamics, which repro- 
duce networks with a selected number of links supporting dy- 
namic criticality by dynamics-correlated addition and deletion 
of links B44 14511 . Similarly to these models, if one allows link 
rewiring only, preserving the total number of links, and make 
it happen depending on the network's dynamical state, the net- 
work is expected to be organized so as to have a broad degree 
distribution, which is under investigation. 

While we focussed on the Hamming distance to capture 
the effects of structural features on the dynamic stability of 
Boolean networks, it would be also interesting to see how 
the heterogeneous connectivity patterns affect the properties 



of attractors in the configuration space, given the recent inter- 
ests 03 El El EH in the scaling behavior of the number 
and length of the attractors in the critical Boolean networks. 
It has been shown that the median attractor length of SF net- 
works is larger than that of networks with fixed in-degree at 
the critical point 1461 RtIi . but much more properties remain to 
be investigated. 

The asymptotic behaviors of the in-degree distributions of 
real transcriptional regulatory networks of E. coli Jit] and 
yeast ill EH are hard to discern due to finite size effects. In 
both organisms, there are about one hundred regulators (nodes 
with outgoing links), which impose a cut-off in the measured 
in-degree distributions. However, one can find for the net- 
work of yeast that its in-degree distribution is much broader 
than that of the networks generated by randomly rewiring the 
links, and that the Hamming distance of the Boolean dynamics 
is much slower than that in the randomized networks, demon- 
strating the contribution of heterogeneous connectivity pattern 
to maintaining dynamic criticality [18]. 

It should be noted that the real transcriptional regulatory 
networks and other biological networks have much richer 
structural properties than described here and their relation to 
the dynamic criticality of the system is of interest. For exam- 
ple, the correlation of the degrees of neighborin g no des has 
been identified in many real- world networks [48, 49] includ- 
ing the yeast gene regulatory network [|20ll and there are stud- 
ies on the effects of the degree-degree correlation on the struc- 
ture and dynamics of complex networks ll50ll5lll52ll . While it 
has been shown [52] that a negative (positive) degree-degree 
correlation is irrelevant (relevant) to the percolation transi- 
tion in complex networks, it still remains to be addressed how 
the correlation affects the critical phenomena of Boolean net- 
works. 
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